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ABSTRACT 

We study the impact of setting initial conditions in numerical simulations us- 
ing the standard procedure based on the Zel'dovich approximation (ZA). As it is 
well known from perturbation theory, ZA initial conditions have incorrect second and 
higher-order growth and therefore excite long-lived transients in the evolution of the 
statistical properties of density and velocity fields. We also study the improvement 
brought by using more accurate initial conditions based on second-order Lagrangian 
perturbation theory (2LPT). We show that 2LPT initial conditions reduce transients 
significantly and thus are much more appropriate for numerical simulations devoted 
to precision cosmology. Using controlled numerical experiments with ZA and 2LPT 
initial conditions we show that simulations started at redshift z% — 49 using the ZA un- 
derestimate the power spectrum in the nonlinear regime by about 2, 4, 8% at z — 0, 1, 3 
respectively, whereas the mass function of dark matter halos is underestimated by 5% 
at m = 10 15 M Q /h (z = 0) and 10% at m = 2 x 10 u M Q /h (z = 1). The clustering 
of halos is also affected to the few percent level at z = 0. These systematics effects 
are typically larger than statistical uncertainties in recent mass function and power 
spectrum fitting formulae extracted from numerical simulations. At large scales, the 
measured transients in higher-order correlations can be understood from first principle 
calculations based on perturbation theory. 

Key words: large-scale structure of the universe - methods: numerical. 



1 INTRODUCTION 

The growth of structure in the universe is a fundamental 
aspect of cosmology, and can be used to constrain the na- 
ture of primordial fluctuations and the energy contents of 
the universe (e.g. dark matter and dark energy). Numerical 
simulations are a very important tool in making theoreti- 
cal predictions of cosmological observables at small scales 
when fluctuations are nonlinear. Fluctuations are followed 
from their primordial seeds through the radiation era to 
the matter dominated era and current acce lerating phase 
by us i ng linear Boltzmann equation solvers ([Peebles fc YrJ 
I l970l: iBond fc Efstathioul Il984t iMa fc Ber^hirrgeiTlF995r 
Isebak&^ZakTirriE^^ which are 

then coupled to N-body simulations to follow their nonlinear 
evolution. This coupling takes the form of final conditions 
predicted from linear Boltzmann codes input as initial con- 
ditions into N-body simulations at some redshift Zi. 

There are competing requirements on the best choice of 
Zi. One would like to start a simulation as early as possi- 
ble, where perturbations are closest to their linear value but 
this can be challenging from numerical reasons since one is 



trying to impose a small perturbation spectrum that can be 
overwhelmed by sources of noise such as numerical round- 
off error or shot noise of the particles. Also, when running 
simulations that include a baryonic component, the initial 
conditions simplify considerably the lower Zi is, in which 
case the baryons had more time to fall into the potential 
wells of the dark matter. 

Despite the rapid progress in the development of differ- 
ent numerical N-body techniques for solving the nonlinear 
growt h of structure in the universe (see e.g. iBertschingeJ 
(1998) for a review), the way initial conditions are input 
into these codes has essentially been the same since the first 
deve lopments in cosmological simulations over twenty y ears 
ago llKlvpin fc Shandarinl | l983t lEfstathiou et alJll985f) . by 
usingthe|ZeFd£vjc hi dl9T0h approximation (hereafter ZA). 

As observations of the late universe quickly approach 
percent-level accuracy, it is vital to examine to what ex- 
tent theoretical predictions can accurately go from primor- 
dial fluctuations to different stages of nonlinear evolution in 
recent epochs, e.g. redshifts z < 3. Studies of linear Boltz- 
mann solvers (ISeliak et al]l2003ft have shown that indepen- 
dent codes agree to within 0.1%, whereas for N-body codes 
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the agreement, while much more difficult to quantify, is of 
the order of few to te n percent depending on the statistics 
and scales considered (|Splinter et alllOQSUFrenk et alJl99i: 
iHeitmann et al.ll2005l) . 

The link between Boltzmann and N-body codes, or how 
to set initial conditions for nonlinear evolution, has been sur- 
prisingly little studied. The accuracy of the ZA has been es- 
sentially taken for granted and instead studies of initial con- 
ditions have concentrated on other (important) issues, re- 
lated to how to best implement the ZA in different situations 
jBertschingej Il995l: ISalmonl Il996t IPenl Il997t iBertschingeJ 
l200ll:lYoshida et alJl2003l : ISirkoll2005lP 

The rationale for this state of affairs seems to be that 
starting a simulation at "high enough" redshift with linear 
perturbation theory (which the ZA is, in Lagrangian space) 
should be safe. A rule of thumb for the choice of starting 
redshift often used is that the rms density fluctuations at 
the interparticle distance scale be of order 0.2. As a result, 
low resolution simulations are typically started much later 
than high resolution ones. 

However, one of the lessons of pertu rb ation the- 
ory (hereafter PT, see IBernardeau et all J200ll) for 
a review) is that the ZA has a very limited va- 
lidity, reproducing only the correct linear growth 
and under estimating the s ke wness and higher-order 
moments jGrinstein fc Wise! Il987t lj uszkie wicz et alJ 
19931: IBernardeau! Il994t ICatelan fc Moscardinil 11994 
Juszkiewicz et alJ Il995t iBaueh et al.1 Il995l : IColombi et alJ 
^9r^) ~with deviations being particularly severe for velocity 
field perturbations because of the ZA failure to conserve 
momentum. 

As a result of this, starting N-body simulations from 
ZA initial conditions leads to incorrect second and higher- 
order growing modes; this excites nonlinear decaying modes 
(called transients) that describe how long it takes for the 
correct dynamics to recover t he actual statistica l properties 
of density and velocity fields iScoccimarrolll998T) . This tran- 
sient evolution is well understood at large scales using PT 
and describes the anomalous behavior seen in measurements 
of higher-order moments in N-body simulations s t arted from 
ZA initial conditions (as shown in IScoccimarrol il998l) and 
demonstrated in much more detail below). 

The importance of transients is not just of interest for 
studies involving higher-order moments, as we show in this 
paper. Underestimation of non-Gaussianity means that the 
tails of the probability distribution of densities and veloci- 
ties will be artificially suppressed, which implies that rare 
events will be even more rare in a simulation that has tran- 
sients. Since typically large volumes (and thus typically low 
resolution) are used to study rare events, the dark matter 
halo mass function at the high-mass end will be particularly 
affected by transients given that low-resolution simulations 
are often started later than high-resolution ones. In addition, 
the transient behavior of large-scale higher-order moments 
implies that nonlinear couplings take a significant time to 
develop their full strength. These are the same couplings 
that act at small scales and give rise to the nonlinear evolu- 
tion of the power spectrum, therefore transients should also 
affect the power spectrum at nonlinear scales. 

In this paper we revisit the issue of transients from 
ZA initial conditions and study numerically the effects on 
higher-order moments, the mass function of dark matter ha- 



los an d the power spectrum. As proposed in IScoccimarrol 
(1998), transients can be reduced significantly by using 
second-order Lagrangian PT (hereafter 2LPT). This leads 
to correct second-order growth factors and greatly improves 
the behavior of higher-order moments. Furthermore, tran- 
sients in this case decay as inverse square power of the scale 
factor, much faster than the inverse of the scale factor in the 
ZA case. We study the improvement brought by using 2LPT 
instead of ZA in the behavior of transients and quantify their 
magnitude in different statistics commonly used. 

The paper is organized as follows. In section [5] we re- 
view how transients arise and present predictions that are 
contrasted with simulations later. Section[3]discusses the nu- 
merical simulations we run to quantify the transients in the 
nonlinear regime. In section |1] we present the main results 
of the paper on the magnitude of transients in higher-order 
cumulants, power spectrum, bispectrum, mass function and 
bias of dark matter halos. We summarize the results and 
conclude in section |^| In appendix ^ we discuss second- 
order Lagrangian perturbation theory and in appendix [B] 
we derive analytically the evolution of perturbations toward 
a glass configuration and discuss glass initial conditions from 
the point of view of transients. 



2 TRANSIENTS 

In this section we discuss why transients exist and how they 
affect the statistics of density and velocity fields at l arge 
scales, following the treatment by IScoccimarrol (|l998h . to 
which we refer the reader for more details. The analytic re- 
sults summarized here also give a framework to understand 
the measurements in simulations we discuss below at small 
scales. Throughout this paper we assume Gaussian primor- 
dial fluctuations. For the most part we assume a flat universe 
with Q m = 1 and Qa — 0, we mention the more general case 
at the end of section 12.11 

2.1 Equations of Motion 

The equations of motion for the evolution of dark matter 
can be written in a compact way by introducing the two- 
component "vector" 

*«(k,f/) = (<5(k, V ), -e(k, n )/n), (i) 

where the index a — 1,2 selects the density or velocity com- 
ponents, with <5(k) being the Fourier transform of the density 
contrast <5(x,r) = p(x)/p — 1 and similarly for the peculiar 
velocity divergence 9 = V- v. Tt = dlna/dr is the conformal 
expansion rate with a(r) the cosmological scale factor and 
r the conformal time. The time variable rj is defined from 
the scale factor by 

n = In a(r), (2) 

corresponding to the number of e-folds of expansion. The 
equations of motion in Fourier space can then be written 
as (we henceforth use the convention that repeated Fourier 
arguments are integrated over) 

ft,*„(k) +n ab *b (k) = 7 $ c (k,ki,k 2 ) *i(ki) *«(ka), (3) 
where 
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a, 



-1 
-3/2 1/2 



(4) 



describes the linear structure of the equations of motion. 
In particular, the first row describes the linearized continu- 
ity equation d^i — ^2 = and the second row describes 
linearized momentum conservation, with the two entries be- 
ing the Poisson equation (V 2 <£> = (3/2)H 2 S) and the Hub- 
ble drag term, respectively. The symmetrized vertex matrix 
~W 

' abc 

Fourier modes and is given by 



^ describes the non linear interactions between different 



7~2 1 (k, ki,k 2 ) = <5 D (k-ki-k 2 
7 222 (k, ki,k 2 ) = <5 D (k-ki-k 2 



(ki + kg) ■ ki 
|ki+k 2 l 2 (ki-ka) 



(5) 



(6) 

7^ ) c (k, ki,k;,) = 7^ s c ) i) (k,k 3 ,ki), and 7 is zero otherwise, S D 
denotes the Dirac delta distribut ion. The formal in tegral 
solution to Eq. is given by f see IScoccimarr J {2000) for a 
detailed derivation) 



* a (k,?7) = g a b(v) <M k ) +/ d v' 9ab(v~v') 



x 7 S(k,k 1 ,k 2 )* c (k 1 ,77')* d (k2 1 r?') I 



(7) 



where a (k) = *l/ a (k,77 = 0) denotes the initial conditions, 
set when the linear growth factor a(r) = 1 and 77 = 0. The 
linear propagator g a b{ri) is simply related to Q a b, and given 
by an inverse Laplace transform (n ^ 0) 



g a b(v) 



ds 
2iri 



(0. ab + s5 ab ) 1 



which given Eq. £Q yields 



g a b(v) 



-377/2 



(8) 



(9) 



This describes the standard linear evolution of the density 
and velocity fields from their initial conditions specified by 
<Mk). 

In the ZA, the Po isson equation is replaced 
by V 2 $ = -(3/2)H 6 ijMunshi fc Starobinskvl Il994t 
iHui fc Bertschingeitil99riF) . and thus Eq. (|3J is replaced by, 



which from Eq. @ gives the linear propagator 



g a b (v) = e " 






1 


+ 


1 


1— 





1 









(10) 



(11) 



The change in the linear propagator (a reflection of gravity 
being driven by the velocity divergence instead of density 
fluctuations) has important consequences, as its structure 
determines the generation of non-Gaussianity. 

From Eq. (J7J we can see why transients arise. For sim- 
plicity, let's assume first that the equations of motion are 
linear, i.e. that j a b c = and the second term in the right 



hand side of Eq. J7J is zero. Then, selecting growing mode 
initial conditions is very simple. All we have to do is to select 
<f} a (k) = <5o(k) u a , where <5o(k) is a Gaussian random field 
whose power spectrum follows from a Boltzmann solver and 
u a — (1,1) selects the growing mode upon contraction with 
g a t in Eq. @. Selecting the growing mode is easy since we 
can solve the linear evolution analytically (given by the first 
term in Eq. 0), otherwise one would generally excite the 
decaying mode. 

In the nonlinear case, we do not know an exact solu- 
tion (that's the purpose of running a simulation), therefore 
we cannot select initial conditions in the exact (nonlinear) 
growing mode. The best we can do is to try to use PT to 
select the fastest growing mode at each order in PT (2LPT 
achieves this to second-order, and does a good job for higher 
orders as well). Since the initial conditions will not be in the 
exact growing mode, decaying modes will be excited by non- 
linear interactions in the second term of Eq. 0: these are 
what we call transients. This is analogous to what would 
happen in the linear case if a (k) was not the eigenvector 
(1, 1), i.e. the decaying mode in the second term of Eq.© 
would be excited. 

In the nonlinear case the problem of transients is worse 
than one might think based on linear theory because tran- 
sients do not go away as fast as in linear theory, where the 
decaying mode is suppresed by a~ 5//2 compared to the grow- 
ing mode. Indeed, if the ZA is used to set initial conditions 
transients are only suppressed by a -1 , so they take a long 
time to decay. Using better initial conditions, such as 2LPT, 
leads to a~ 2 suppression, a substantial improvement. 

To describe transients in a quantitative way, we expand 
the initial and final conditions in terms of Gaussian linearly 
evolved density perturbations <5o(k), 

00 00 
Q (k)=^^(k), * Q (k, J7 )=^*i" ) (k, ?7 ), (12) 

n— 1 n— 1 

with 

0i n) (k) = [5v]n Xi n) (ki, . ..,k„) J (ki) . ..SoQtn), (13) 

*i n) (M) = [8v] n T a n) ^...^ n ;n) <5o(ki)...Mkn),(14) 

where [<5 D ]n = &>(k — ki...„) and ki...„ = ki + . . . + k„. The 
kernels X a n ' describe how the initial conditions are set (at 
17 = Q), e.g. if linear theory is used = u a = (1, 1) and 
X a "^ — for n > 1. The kernels J- a n \rj) describe how per- 
turbations evolve after initial conditions are set, they obey 
the recursion relation jScoccimarrolll998l) . 

m=l u 



+ g ab (v) 4 n) . 



(15) 



obtained by replacing Eq. iTEZl and Eqs. 1131141 into Eq. J7J. 
For brevity we have suppressed the Fourier arguments. Note 
that J- a n \n = 0) = X„ as it should be, afterwards the 
fastest growing mode behaves as T a (rj) ~ e nv at each or- 
der n. Transients determine how fast one goes from initial 
conditions given by I a n ^ to the regime where final conditions 
are in the fastest growing mode Ta (n) ~ e n7) . 
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If the initial conditions are not set using linear theory, 
the kernels T a obey recursion relations that follow from the 
dynamics used to set them. For ZA initial conditions T a 
obeys, 



r(«) 



n-l 

E 

m— 1 



(16) 



with the condition li 1 ' = (1,1). This can be understood 
from Eq. I15H by using linear theory initial conditions and 
keeping the fastest growing mode, formally done by evolv- 
ing by an infinite amount of time (from s = — oo). Note that 
these kernels are independent of rj, that's guaranteed due to 
the lower limit in the time integral. In the long time limit, 



the T a (rj) kernels approach Ta (r)) = e" v H a , where 
Ha obeys the recursion relation in Eq. 1161 with the lin- 
ear propagator replaced by g ao , Eq. It is easy to 
check that if Zl™' is replaced by such Ha in Eq. (1151 , then 



is a solution at all times, i.e. there are no 
transients, as expected since one is setting initial conditions 
exactly, i.e. in the fastest growing mode. 

Although most simulations are started by initial con- 
ditions decribed b y the kernels in Eq. (1 1 6} when using the 
ZA, some authors lEfetathiou et alll98alPower et 811 2003) 
define the velocities in a different way, by applying the ZA 
prescription to the perturbed density field after particles have 
been moved inst ead of the linear d ensity field. In this case, 
it can be shown dScoccimarrolll99sl) that the velocity diver- 
gence kernels are not those in the ZA, but rather they are 
equal to the density kernels in the ZA. That leads to a re- 
duct ion in the time sca le of transients by a factor of about 
two iScoccimarrdll998l) . We will not consider this possibil- 
ity here, as the standard ZA starts are most often used in 
the literature, and the improvement brought by using this 
prescription for the velocities does not fully address the lim- 
itations of the ZA (e.g. the second-order growth is incorrect, 
and the transients still evolve as the inverse of the scale fac- 
tor). 

As an example, the second-order density kernel reads 
from Eq. 1151 , with x = k\ • k^, 



/(") 



+ c 



+ 



!( 



ti _|_ hi 

&2 ki 



2 2 



77 3 , 2 

e V* 



-3V/2 J, x 2 

35 v 



which interpolates between the ZA value at rj = and the 
exact value for large ij (that in between square brackets). 
Note that, as mentioned above, the dominant transient term 
is suppressed only by a single inverse power of the scale fac- 
tor, so it takes a long time to disappear. If 2LPT is used 
to set initial conditions, the second-order kernels are repro- 
duced exactly and there are no transients to this order. For 
higher-order kernels, 2LPT shows trans ients suppressed b y 
inverse square powers of the scale factor dScoccimarrcl l998'). 
and as we see below they are greatly reduced compared to 
the ZA case. 

In making predictions that evolve all the way to z = 
one must take into account the dependence on cosmological 
parameters neglected so far. A simple and rather accurate 



way of doing this is to redefine in Eq. Q, —6(li,r))/fH as 
the second component, where / = d In D+ jd In a and D+ is 
the growth factor with the exact cosmological dependence. 
When this is done, the only change in the equations of mo- 
tion is 3/2 — » 3f2 m /2 / 2 in Eq. (|1J, and given that / ~ Qm 9 
jBouchet et alJll995h for flat models with a cosmological 
constant this can be well approximated by the TO = 1 value 
of 3/2 since during most of the time evolution Q m /f 2 is 
very close to unity. This leads to the same kernels as in the 
fl m = 1 case but with the scale factor replaced everywhere 
by the growth factor D+. 



2.2 Statistics 

We now discuss quantitatively the impact of transients on 
statistics of the density field at large scales. Since linear 
growing modes are reproduced exactly, the power spectrum 
at large scales is not affected by transients. The effect of 
transients a large scales is only felt by higher-order spec- 
tra, such as the bispectrum, or cumulants, e.g. the skew- 
ness and kurtosis. However, as these quantities characterize 
non-linear couplings, they are useful for understanding the 
results at small scales where nonlinearities play a role in the 
power spectrum, which will then be affected by transients, 
as we shall see using simulations in section [4~2*1 

The power spectrum -P(fc), bispectrum B123 = 
B(ki, fc 2 , ^3) and trispectrum T1234 are given by, 



<*(k)*oo> 

(*(ki)*(ka)*(ks)) 
<5(k 1 )5(k2)<5(k 3 )<5(k 4 )) c 



<5 D (k + k') P(k), 

<5D(kl23) -Bl23, 
<5D(ki234) 7l234, 



(18) 
(19) 
(20) 



with ki23 = ki +k2 + k3, etc., and the subscript "c" denotes 
a connected contribution, whereas the skewness 5*3 and kur- 
tosis S4 are given by, 



S3 = 



(s 3 ) 



S 4 



(21) 



In second-order perturbation theory the bispectrum reads 



B123 = 2JFf ) (ki,k 2 ;T?) P (fci)Po(fc2) + cyclic, 



(22) 



where Po denotes the power spectrum of the linear fluctua- 
tions at the initial conditions (i.e. P(k) = a 2 P (k) in linear 
theory), and cyclic permutations over {fci, fe, kz\ are under- 
stood. The first three cumulants are given by, 



(17) (S 2 ) 



(S 3 ) 



P(k) W 2 (kR) d A k, 



-Bl23 W1W2W3 fc(ki23 )d 3 k l d 3 k 2 d 3 k 3 , 



(23) 



(24) 



{8 4 }c = / T 123 4 WW2WW4 fo(kias) d 3 fci . . . d 3 k 4 (25) 



where W(kR) represents a top hat filter in Fourier space at 
smoothing scale R, and W% = W(kiR). From these expres- 
sions i t follows the skew ness and kurtosis including tran- 
sients iScoccimarrd H998'l 



(34 1 l r , , 1/ 26 
S 3 = ( y - 7l } + -[4- 7l ] + -( 7l - T 



12 



35a 7 / 2 ' 
(26) 
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Figure 1. The skcwness S3 at R = 10 h~ 1 Mpc for Zel'dovich 
approximation (solid) and Eulerian linear theory (dashed) initial 
conditions as a function of the expansion factor a away from initial 
conditions (at a = 1). Values are given in terms of the exact value 
of S3. 



S 4 = 



60712 _ 

1323 
1312 

+ 



6271 7t? 
3 3 
871 



35a 5a 
192 

245a 7 / 2 5a 7 / 2 4725a 9 / 2 + 1225a 7 ' 



272 
3 

1504 



816 2871 



+ 



+ 



184 
75a 2 

(27) 



where we have written the time dependence in terms of 
the scale factor a and 



7p 



d p lna 2 (R) 
d\n v R ' 



(28) 



with a 2 (R) the variance of density perturbations, Eq. 1231 . 

These results have also been derived more recently using 
a different approach based on the ste epest descent method 
for the density PDF JValageasll2002T) . In Eq. we have 
explicitly written down the different contributions: the ZA 
contribution from initial conditions in square brackets, the 
exact value (reached only asymptotically as a — > 00) in curly 
brackets, and the remaining two transient terms (which can- 
cel the exact value when a = 1). If initial conditions were 
set using Eulerian linear theory, the term in square brackets 
would not be present. 

As discussed above, to include the dependence of tran- 
sients on cosmological parameters we will make the replace- 
ment a — > D+ in all comparisons of Eqs. 126I27H with nu- 
merical simulatons. 

FigureQevaluates Eq. I26H in two cases, one where ZA is 
used to set the initial conditions (solid) , and the other when 
Eulerian linear theory is used (and thus S3 = at a = 1). We 
see that if linear theory were used, to reach the correct value 
of S3 within some tolerance, it would take about four times 
more expansion away from the initial conditions compared 
to the ZA case. The difference in the two cases is purely 
due to the initial value of S3, i.e. the rate of convergence to 
the right answer is the same (as transients go away as a -1 



in both cases). In contrast, using 2LPT initial conditions 
improves the initial values of the S p parameters dramatically 
over the ZA case (being exact for S3), and also improves 
substantially the rate of decay of transients to aT 2 . 

Linear theory initial conditions are relevant for codes 
that do no use particles, e.g. in the case where baryons are 
described by fields in a grid instead of particles. In such cases 
one must keep in mind that the magnitude of the transients 
we demonstrate with simulations below are very likely to 
be much worse, as reflected by the difference between the 
curves in Fig. 

The dependence of transients on the shape of the power 
spectrum, in Eqs. 12612711 through y p , shows that transients 
from ZA initial conditions are somewhat reduced as the 
spectrum becomes steep, since the gravitational interactions 
are dominated by large scales and this leads to fewer deflec- 
tions, i.e. more coherent displacement fields better approxi- 
mated by straight-line trajectories given by the ZA. The de- 
pende nce on spectrum is discussed in detail in IScoccimarrol 
(1998). Here we concentrate on a fixed shape given by a 
concordance cosmological model as discussed in the next 
section. 

The lesson from Fig. is that nonlinear couplings (in 
this case, second order) take a significant time to develop 
their full strength if initial conditions are not set accurately. 
As expected from perturbation theory, these weaker nonlin- 
ear couplings give a transient effect in the nonlinear power 
spectrum. The effect of transients is also felt rather strongly 
on rare events such as the abundance of high-mass clusters, 
which probe the tail of the density probability distribution 
characterized by the S p parameters. We demonstrate these 
effects below using numerical simulations. 



3 SIMULATIONS 

In order to study the effects of transients at small scales we 
ran a set of simulations with different initial conditions (ZA 
and 2LPT), starting redshifts z%, box sizes Lbox, softening 
lengths and number of timesteps as described in Tabled All 
our simulations have iVp ar = 512 3 particles with cosmologi- 
cal parameters Q, m = 0.27, fi A = 0.73, Q b = 0.046, h = 0.72 
and as(z = 0) = 0.9, and were run using the Gadget2 code 
jSpringell2005l) . 

All the internal parameters to the Gadget2 code used to 
run the different numerical experiments are also given in Ta- 
blc0 In summary, we have two sets of initial conditions (ZA 
and 2LPT), each of them started at two different redshifts 
Zi with the purpose of checking the magnitude of transients. 
The difference between ZA and 2LPT initial conditions is 
explained in Appendix^ in terms of Lagrangian displace- 
ments. In our implementation, particles are started from a 
grid and all differentiations are done in Fourier space, i.e. 
through V — > ik. Another possibility would be to start par- 
ticles from a "glass" , in appendix[B]we briefly discuss how to 
avoid inducing extra initial skewness (and thus transients) 
due to glass initial conditions. 

A simulation without transients should have statistical 
properties (measured at redshift z < «,) independent of z%. 
Also, if generating initial conditions with the ZA was ac- 
curate enough, certainly 2LPT would be as well, since PT 
converges very well at high redshift. If this were the case, no 
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Table 1. Numerical simulations parameters. All our simulations have N paT = 512 3 particles with Q m = 0.27, 
n A = 0.73, U b = 0.046, h = 0.72 and <r 8 (z = 0) = 0.9, were run using Gadget2 iSpringeJ I2OO5I) with time 
integration accuracy parameters ErrTolIntAccuracy=0.025, MaxRMSDisplacementFac=0.2, MaxSizeTimestep=0.025, 
and gravitational force accuracy parameters ErrTolTheta=0.5, TypeQf 0peningCriterion=l, ErrTolForceAcc=0.005 
except for the HQ runs, where more stringent parameters (ErrTolIntAccuracy=0.01, MaxRMSDisplacementFac=0.1, 
MaxSizeTimestep=0.01, ErrTolTheta=0.2, ErrTolForceAcc=0.002) were chosen. Length units are in h~ 1 Mpc and 
mass units are in 10 10 Mq//i. 



Name 


^box 




Initial Conditions 


Zi 


-^realizations 


softening 


timesteps 


LRZA49 


1024 


59.94 


ZA 


49 


8 


0.2 


800 


LRZA24 


1024 


59.94 


ZA 


24 


8 


0.2 


750 


LR2LPT49 a 


1024 


59.94 


2LPT 


49 


8 


0.2 


800 


LR2LPT11.5 


1024 


59.94 


2LPT 


11.5 


8 


0.2 


650 


LRsZA49 


1024 


59.94 


ZA 


49 


2 


0.07 


1800 


LRs2LPT49 


1024 


59.94 


2LPT 


19 


2 


0.07 


1800 


HRZA49 


512 


7.49 


ZA 


49 


4 


0.04 


2200 


HRZA24 


512 


7.49 


ZA 


24 


4 


0.04 


2200 


HR2LPT49 6 


512 


7.49 


2LPT 


49 


20 


0.04 


2200 


HR2LPT11.5 


512 


7.49 


2LPT 


11.5 


4 


0.04 


2200 


HQZA49 


512 


7.49 


ZA 


49 


1 


0.02 


6500 


HQ2LPT49 


512 


7.49 


2LPT 


49 


1 


0.02 


6500 



a This is the reference run for Lb ox = 1024. 
^ This is the reference run for Lbox = 512. 

appreciable difference would be seen at low redshift between 
ZA and 2LPT initial conditions. 

However, for ZA initial conditions, we observe that 
statistics at low redshift are different than those from 2LPT 
initial conditions started at the same redshift (in this case 
Zi = 49), and we also see a clear dependence on Zi (com- 
paring Zi — 49 with Zi — 24 simulations). Furthermore, 
the magnitude of transients observed at large scales agrees 
with the perturbative predictions presented in the previous 
section. In addition, the expected suppression of nonlinear 
couplings at large scales acts at small scales in a manner 
consistent with what one expects physically. 

How do we know the 2LPT initial conditions simula- 
tions are correct? At large scales we checked the numerical 
simulations against the PT predictions, and we see almost 
no dependence of the low redshift results on the starting red- 
shift Zi. To make this latter point more obvious, we decided 
to present results for a very late start, Zi = 11.5 (compared 
to our standard choice Zi = 49). Note we are not advocat- 
ing starting simulations at Zi = 11.5, this is only to prove 
that 2LPT initial conditions have very little transient effects. 
2LPT initial conditions do have transient effects, although 
much smaller than ZA. Starting at redshift z% = 11.5 makes 
them large enough to be measurable at z — 3, which allows 
us to estimate the magnitude of transients in our reference 
runs (see section |K| for a discussion of this). 

Most of the results presented here are obtained by using 
the "high-resolution" runs (labeled HR in Table , corre- 
sponding to Lbox = 512fo _1 Mpc. The "low-resolution" (la- 
beled LR) simulations with Lbox = 1024 ft -1 Mpc are used 
to study the high-mass tail of the dark matter halo mass 
function in Section [4.51 In addition, we have studied if the 
magnitude of transients depends on the softening length by 
reducing the softening length by a factor of three in the LR 
runs (denoted by LRs), and we have also run HR simula- 



tions with twice smaller softening length and stricter time- 
integration and force accuracy parameters, labeled "high- 
quality" (HQ) runs in Table Q From these sanity checks 
we conclude that the magnitude of transients observed from 
ZA initial conditions we present below is insensitive to the 
choices we have made in the less time-consuming runs. 

The runs with 2LPT initial conditions with Zi = 49 
are our most accurate simulations from the point of view of 
transients. We denote such realizations as reference runs. In 
order to quantify the magnitude of transients, we measure 
statistics at low redshift and compare them to the same 
measurements in these reference runs. 



4 RESULTS 

4.1 Higher-Order Cumulants 

We first start by comparing simulations for measurements 
of higher-order cumulants at large scales, skewness 53 and 
kurtosis Si, which are sensitive to transients and can be 
understood analytically using perturbation theory, as pre- 
sented in Section f2. 21 

Figure[2]shows the measurement of S3 and S4 as a func- 
tion of smoothing scale R for the reference runs HR2LPT49, 
after averaging over 20 realizations. The agreement at large 
scales with the perturbative prediction (solid lines) is ex- 
cellent, better than 1% for S 3 for R > 40 h' 1 Mpc. Note 
that all measurements at different redshifts z = 0,1,2,3 
lie on top of eac h other in the large-scale limit as predicte d 
by tree- level PT l|.Tuszkiewicz et, al.ll993l : lBernardeatJl994T ). 
and depart from it at increasingly l arge scales as z decreases 
as ex pected from one-loop PT l|Scoccimarro fc Friemanl 
Il99fit iFosalba fc Gaztanaeal ll99fiT l . Previous measurements 
in large-volume simulations such as the Hubble volume 
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Figure 2. The skewness S3 (bottom) and kurtosis S4 (top) for 
the reference runs (2LPT Zi = 49) at redshifts 2 = 0,1,2,3 (cir- 
cles, pentagons, squares, and triangles, respectively, from top to 
bottom in each panel. 
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Figure 4. Same as Fig. 151 but for ZA Zi = 24 initial conditions. 
The solid lines show the predictions of the expected transient 
behavior in S3, Eq. and S4, Eq. 1271 . 
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Figure 3. The skewness S3 (bottom) and kurtosis 54 (top) for the 
2LPT Zi = 11.5 initial conditions compared to the reference runs 
(zi = 49) at redshifts z = 0, 1, 2, 3 (circles, pentagons, squares, 
and triangles, respectively, from top to bottom in each panel). 



iColombi et al.l2000l) could not achieve such a clean measure 
of higher-order cumulants independent of z due to transients 
effects from Zi — 35 ZA initial conditions, as will become 
clear shortly. 

Figure 13 shows what happens with 2LPT initial condi- 
tions when the starting redshift is changed to a very late 
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Figure 5. Same as Fig. l3lbut for ZA Zi = 49 initial conditions. 
The solid lines show the predictions of the expected transient 
behavior in S3, Eq. (2<3, and S4, Eq. l27l . 



start, Zi — 11.5. As expected theoretically, at large scales 
there should be no change whatsoever in 53 since 2LPT 
reproduces exactly the second-order growing modes, and 
that's indeed what is seen in the bottom panel. At third- 
order, 2LPT does underestimate 5*4 slightly, and that mani- 
fests in a small transient behavior at large scales (top panel). 
Note that since third- and higher-order couplings are slightly 
smaller in 2LPT, one expects the loop corrections to 53 and 
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Figure 6. Power spectrum for different initial conditions (2LPT 
Zi = 11.5, ZA Zi = 49 and ZA Zi = 24, from top to bottom in 
each panel) compared to the reference runs at z = (top), 2 = 1 
(middle) and z = 3 (bottom). The dotted lines show an estimate 
of the transients for ZA Zj = 49 and ZA Zi = 24 from one-loop 
PT. 



Figure 7. Rcdshift-space power spectrum for ZA Zj = 49 initial 
conditions compared to the reference runs at z = (top), z = 
1 (middle) and z = 3 (bottom). In each panel we show modes 
parallel to the line of sight (triangles), monopole (circles) and 
modes perpendicular to the line of sight (squares). 



S4 to be smaller than in the reference runs, that's indeed 
what is seen in both panels in Fig. El at small scales: there 
is a slight underestimate of the small-scale cumulants, al- 
though by the time one reaches 2 = 0, that is well under 
1%. Overall, the behavior of 2LPT initial conditions is re- 
markably stable given the very late start, Zi = 11.5. In fact, 
we picked such a late start to make the effects visible at all 
at 2 < 3. 

Figures 01 and El show the corresponding situation with 
ZA initial conditions for Zi = 24, 49, respectively. It is obvi- 
ous from these figures that the ZA initial conditions runs are 
not stable at large scales, in particular, z% = 24 and Zi = 49 
ZA runs do not give the same results, with significant de- 
viations at the few percent level. These deviations are in 
excellent agreement with the PT predictions, Eq. 1261 for 
the skewness (bottom panels) and Eq. 1271 for the kurto- 
sis (top panels) . From this we conclude that the simulations 
started with ZA initial conditions have transients which are 
well understood at large scales from first principles. Figs. El 
El check the PT transients pred ictions more accurately than 
done before <IScoccimarro]ll99sl) . We now discuss how these 
effects impact on clustering at small scales. 

4.2 Power Spectrum 

The power spectrum is the most widely used statistic. Non- 
linear corrections to the power spectrum are controlled by 
the same nonlinear couplings that determine higher-order 
statistics at large scales, and thus we expect to see transient 
behavior in the nonlinear power spectrum as well. 

Figure|S|shows that this is indeed the case. First, in each 
panel the top line denotes the 2LPT 2; = 11.5 initial condi- 



tions runs, showing that in the 2LPT case there is almost no 
transient behavior, as expected from the results on the S p 
parameters at large scales. On the other hand, the ZA initial 
conditions runs show significant transient effects, up to 10% 
at z = 3, which can potentially lead to serious systematic 
errors for precision studies of the high-redshift universe. For 
example, a ZA start at z< = 30 was used for the LCDM 
simulation from which the latest fitting fo rmula for the non- 
linear power spectrum ijSmith et alJl2003h was derived. Our 
reference runs show a power spectrum for k < 2/iMpc -1 
that is as much as 7% higher at 2 = and as much as 
13% h igher at 2 = 3 when compared to the ISmith et alJ 
(2003) fitting formula. These deviations are in reasonable 
agreement with the expectations based on Fig. El A detailed 
comparison of the nonlinear power spectrum in our simula- 
tions to fitting formulae, one- loop PT and renormalized PT 
flCrocce fc ScoccimarrolfeoOfit ) will be discussed elsewhere. 

It is also important to note that a relatively high- 
redshift ZA start such as z% = 149, which we ran as a test, 
leads only to mild improvements compared to ZA Zi = 49 in 
Fig. El i.e. the suppression in the nonlinear regime is still of 
order 1%, consistent with the a -1 slow scaling of transients 
in the ZA. We have checked that these results do not depend 
on the accuracy of the force, time integration, and soften- 
ing. For example, using the HQ runs (see tabled for 2; = 49 
we obtain very similar results to those in Fig. El specifically 
within about 0.04% for z = 0, 1 and 0.02% for 2 = 3. Also, 
as shown in Fig. El the magnitude of the measured transients 
in the ZA z% = 49 and ZA Zi = 24 runs are in reasonable 
agreement with their estimation by one-loop PT using the 
kernels in Eq. 1)15(1 . particularly at high redshift where the 
agreement should be best. 

Figure |7| presents results in redshift space for the ZA 
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Figure 8. Bispectrum for different initial conditions (2LPT Zj = 
11.5, ZA Zi = 49 and ZA Zi = 24, from top to bottom) compared 
to the reference runs at z = 0. The triangles shown correspond 
to fci = 0.12 /iMpc -1 , &2 = 2fci with 9 the angle between ki and 
k 2 . 



Zi — 49 initial conditions runs. The combined effect of den- 
sity and velocity trans i ents l eads to nontrivial behavior. As 
shown bv IScoccimarrd il998l) transients in the velocity field 
are larger than in the density field. For modes parallel to 
the line of sight (triangles in each of the panels in Fig. 0, 
small-scale velocities suppress power due to velocity disper- 
sion. Therefore, although transients in real space make the 
power smaller, when mapped into redshift-space with veloc- 
ities that are even more suppressed by transients than the 
density, the overall effect is to increase the power spectrum 
along the line of sight compared to the reference runs that 
have no transients. Since there is cancellation between den- 
sity and velocity transients, the overall effect is somewhat 
weaker than in the real space power. 

For modes perpendicular to the line of sight, on the 
other hand, the effect is opposite since these modes are 
unaffected by redshift-space distortions and show the same 
suppression in power seen in Fig. HJ This means that when 
the redshift space power spectrum is averaged over angles 
with respect to the line of sight to obtain the monopole, 
these opposite behaviors lead to cancellation that makes the 
monopole less sensitive to transients (see circles in Fig. 0. 
However, the anisotropy of the redshift-space power spec- 
trum is affected by transients much more than the monopole. 



4.3 Bispectrum 

The bispectrum is also affected by transients, and from 
Eqs. JT7J and we see that at the largest scales where 
second-order Eulerian PT holds it should be affected most 
for triangles that are most different from elongated triangles 
(which correspond to x = 1 and show no transient behav- 
ior) . Figure |H| shows the bispectrum at z = for different 



Figure 9. The PDF of the density contrast <5 in units of rms 
value, a = {5 2 ) 1 / 2 , for our reference run (2LPT Zi = 49) and ZA 
Zi = 24 initial conditions. 



initial conditions divided by that in the reference runs for 
triangles with ki — 0.12 h Mpc -1 and k% = 2k\ as a func- 
tion of the angle 9 between ki and k.2. From this we see 
again that the 2LPT z% = 11.5 initial conditions show very 
little (< 0.2%) transient effects, whereas the ZA initial con- 
ditions runs show the expected transient signature, minimal 
at elongated triangles and maximal at isosceles triangles. 
These can lead to systematic errors in the determination of 
bias and cosmologic al parameters, given present expected 
observational errors JSefusatti et al.ll2005l) . 



4.4 The PDF 

The suppression of higher-order moments by transients 
means that the probability distribution function (PDF) of 
density fluctuations is also affected. Figure [5] shows a com- 
parison of the PDF of our reference run (2LPT z% = 49) and 
ZA Zi = 24 initial conditions at redshift z — 3 for a smooth- 
ing scale R — 10.84 ft -1 Mpc. Note how the non-Gaussian 
features of the PDF are changed by transients: the right tail 
is suppressed, whereas the left cutoff is enhanced, leading to 
smaller skewness. Figure flUl shows the ratio of the PDF for 
different initial conditions to that of the reference run, to 
better appreciate the differences. Although the differences 
at the left tail are more significant, the PDF is falling very 
steeply, thus in practice the differences seen at the right tail 
are more important. These correspond to high-density re- 
gions, and thus we expect these differences to impact the 
high-mass tail of the mass function of dark matter halos. 

4.5 The Mass Function 

Figures ITTI and WI\ show the ratio of mass functions for dif- 
ferent initial conditions to the reference runs, where the 
LR simulations have been used at the high-mass end to 
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Figure 10. Ratio of PDF's for different initial conditions to our 
reference run. Solid lines denote ZA Zi = 24 at z = 3, dashed 
lines ZA z; = 49 at z = 3, dot-dashed ZA z; = 24 at z = 0, long 
dashed lines ZA z; = 49 at z = 0. 



Figure 12. Mass function for ZA initial conditions with z; = 49 
(circles) and z; = 24 (squares) compared to the reference runs at 
z = 1. 
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Figure 11. Mass function for different initial conditions (2LPT 
Zi = 11.5, ZA Zi = 49 and ZA Zi = 24, from top to bottom) 
compared to the reference runs at z = 0. 



probe halos with masses larger than m ~ 10 14 Mq/H. All 
masses have been c orrected by the procedure described in 
IWarren et alJ J2005) to take into account the finite number 
of particles making up halos. We use the friends-of- friends 
algorithm with linking length parameter equal to 0.2 to iden- 
tify halos. 

The effect in the suppression of the high-mass tail is 



clearly seen in the simulations that were started with ZA 
initial conditions. The mass function in ou r reference runs 
is 15% larger than the lWarren et all <l2005l) fitting formula 
for m = lQ^h^My at z = (a nd thus about 30% larger 
than the ISheth fc Tormenl ll 19991) mass func tion), whereas 
for m < 10 14 fe ~ 1 M ( 7) is consi s tent with the IWarren et alJ 
J2005tl fit. The IJenkins et all (1200 if) fit is very similar to 
the I Warren et aTT(l2005T) fit in the mass range covered by 
our simulations. These differences, in good agreement with 
the deviations seen in Fig. 1111 are thus most likely due to 
transients in the simulations used to derived these fitting 
formulae at the high-mass end, which were started from Zj = 
24, 35 ZA initial conditions. 

Such a suppression of the mass function is clearly impor- 
tant for current and future observations that use the abun- 
dance of massive clusters to probe cosmology. It is also worth 
noting that the suppression of the mass function high-mass 
tail due to transients can be made worse by requiring that 
simulations start at a redshift Zi such that density fluctua- 
tions at the interparticle distance be a fixed number (often 
chosen to be about 0.2). This means that the large volume 
simulations that are required to probe rare events will be 
started later and thus will be subject to more suppression 
due to transients. 

Using the LRs simulations with smaller softening length 
(see table for Zi = 49 we obtain very similar results to 
those in Fig. Illl for z — 0, specifically the average deviation 
in the eight bins with m > 10 14 Mq/H (excluding the last 
one that has substantial errors) is within 0.7% of that in 
Fig.HU Therefore the suppression seen in the mass function 
is r obust to a rather l arge change in softening length. 

iReed et alJ ( 2003) studied the effects of the starting red- 
shift Zi on the mass function measured at high redshifts 
z ~ 7 — 15 from numerical simulations initialized with the 
ZA. Their statistical errors are rather large, but they do 



© 2006 RAS, MNRAS OOO.ITllTll 



Transients from Initial Conditions in Cosmological Simulations 11 




2xl0 14 < m < 10 1 
I i L 



8 10 
r [Mpc/h] 



20 



Figure 13. Ratio of halo-halo correlation functions £hh at z = 
for different initial conditions, 2LPT Zj = 11.5 (triangles), ZA 
Zi = 49 (circles) and ZA Zi = 24 (squares), compared to the 
reference runs at z = 0. The different panels correspond to three 
bins in halo mass as indicated (masses are in units of /i -1 Mq). 



seem to see the effects we report here when comparing sim- 
ulations with Zi = 69, 139, 279. They conclude Zi — 139 is 
safe for measurements at z ~ 7 — 15 and find a suppr ession 
of the high-mass tail from the lSheth fc TormerJ Jl999l) mass 
function which becomes stronger with z. Extrapolating our 
results to their range, we would expect an important con- 
tribution from transients for those choices of {z,Zi} going 
exactly in the same direction with increasing z. Therefore 
we suggest that high-redshift studies of the high-mass tail of 
the mass function should be carefully checked for transients. 



4.6 Halo Bias 

Finally, Fie. ll3l shows how the clustering of dark matter ha- 
los is affected by transients at z = 0. The panels show the 
halo-halo correlation function £hh for different initial con- 
ditions compared to those in the reference runs, for three 
different mass bins. Here we have used HR simulations for 
the first bin in mass, and the LR simulations for the two 
largest mass bins. Note that as halo mass becomes larger, 
corresponding to rarer events, the halo bias in the simula- 
tions with more transients is larger. This is consistent with 
the fact that transients suppress the abundance of rare halos 
making them more rare and thus more biased, as expected 
from the peak-background split. 



5 SUMMARY AND CONCLUSIONS 

In this paper we considered the use of different approxima- 
tions to set up initial conditions for cosmological simula- 
tions, concentrating on typical initial redshifts Zi ~ 50, and 
studied their influence on low-redshift (z < 3) statistics. The 



standa rd procedure in the literature is to use the lZerdovichl 
(1970) approximation (ZA), a linear solution in Lagrangian 
space, or even linear perturbation theory in Eulerian space 
(typically when simulating baryons using methods which are 
not based on particles). 

At high redshifts (z ~ 50) the naive expectation is that 
linear perturbation theories should be accurate enough, but 
they are not. The reason for this is that at high redshifts non- 
linear corrections are relatively much more important than 
at low redshift because the relevant scales involved have a 
very steep spectrum (with the effective spectral index very 
close to n c ff = —3). An example of this is shown in Fig. |3] 
where the density PDF smoothed on 10/i _1 Mpc at z = 3 
shows strong non-Gaussian features, even though naively 
these scales are in the linear regime (since the linear rms 
density fluctuation is only a = 0.25). Modeling these non- 
linear corrections using the ZA leads to a significant under- 
estimate of nonlinear couplings and thus non-Gaussianity. 
This leads to a suppression of the density PDF tail for large 
densities that translates into a suppression of the high-mass 
tail of the mass function. The weaker nonlinear couplings 
lead to a smaller power spectrum in the nonlinear regime. 

The mismatch between the approximation used to set 
initial conditions and the correct dynamics leads to tran- 
sients, excitation of decaying modes that are long-lived and 
survive until low redshift, systematically biasing low the 
nonlinear couplings even after evolving with the correct dy- 
namics of an N-body code for as much as 50-100 scale fac- 
tors. Transients cannot be avoided, but they can be signif- 
icantly reduced by imposing initial conditions with better 
approximations to the nonlinear equations of motion. 

More specifically, in this paper we have shown that, 

(i) The measurements of statistics at low redshift for sim- 
ulations started with ZA initial conditions depend signifi- 
cantly on the starting redshift Zi, signaling the presence of 
transients. 

(ii) This can be avoided to a large extent by using 2LPT 
initial conditions, which as we have demonstrated show very 
little dependence on Zi even for extremely late starts such as 
Zi = 11.5. Given that such runs show at most 1% transient 
effects and that transients in 2LPT decay as inverse square 
of the scale factor, we can estimate that our reference runs 
with Zi = 49 are accurate (in terms of transients) to better 
than a tenth of a percent for z < 3. 

(iii) The simplest way to quantify transients is to follow 
how non-Gaussianity develops at large scales as nonlinear 
couplings reach their full strength during the transient evo- 
lution. Our results for higher-o rder cumulants at large scales 
confi rm previous PT results iScoccimarrd Il995 IValaeea3 
120021) with much better accuracy than done before. 

(iv) The small scale behavior of transients is consistent 
with the reduction of nonlinear couplings observed at large 
scales, and we have checked it is not affected by the choice 
of softening length, force and time integration accuracy. 

(v) The use of the ZA at Zi = 49 leads to suppressions 
in the power spectrum of order a few (z = 0) to almost ten 
percent (z = 3) in the nonlinear regime. These results are 
in reasonable agreement with an estimation of transients in 
the power spectrum based on one-loop PT. 

(vi) ZA initial conditions also lead to a suppression in the 
high-mass tail of the mass function of dark matter halos. For 
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Zi = 24 this suppression is about 10% for m — 10 15 h~ Mq 
at z — and m = 10 14 h -1 Mq at z = 1. The clustering of 
halos is also affected at the few-percent level. 

An obvious alternative to using 2LPT initial conditions 
is to start with the ZA much earlier than Zi ~ 50. However, 
in order to compete with the improvement brought by using 
2LPT very high redshifts are required. For example, Fig. [(J 
shows that az, = 11.5 2LPT start has approximately a 0.5% 
transient effect on the nonlinear power spectrum at z — 0, 
whereas a Zi — 49 ZA start has a 2.5% effect. In order for 
the ZA to reduce the transients to the level of the 2LPT 
Zi = 11.5 start, a zi — 249 ZA start would be required. 
Reaching the level of a 2LPT Zi — 49 start, corresponding 
to our reference runs, will require about another factor of 
16 in redshift (given that transients scale as a~ 2 in 2LPT 
compared to a -1 in the ZA). Other issues make starting 
simulations at very high redshifts a bit more challenging, e.g. 
depending on the N-body code algorithm, numerical errors 
may be more difficult to avoid at very high redshifts where 
perturbations are so small. In addition, initial conditions for 
baryons at such high redshifts must be carefully included as 
they are still falling into dark matter potential wells. 

Although we have not studied this in any detail, simu- 
lations started with linear theory instead of the ZA are ex- 
pected to be much more significantly affected by transients, 
e.g. a simulation started with linear theory at z% = 49 should 
be roughly equivalent to one started with the ZA at z = 11.5 
(see Fig. . Note that linear theory initial conditions are 
often used in codes that rely on grid methods with gas dy- 
namics solvers. We expect such methods to show differences 
with e.g. SPH methods that use ZA initial conditons. When 
the dynamical variables are densities and velocities in a mesh 
(instead of information carried by particles), a way to reduce 
transients would be to use higher-order Eulerian perturba- 
tion theory. 

While in this paper we have concentrated on statistics 
of the dark matter density perturbations, it would be inter- 
esting to know how much do transients affect other high- 
redshift probes, where transients are most likely to play 
an important role. In this regard, the most pressing aspect 
would be to check the impact of transients in the Lyman- 
a forest power spectrum, because of the high accuracy de- 
manded from simulatio ns to obtain cosmological constraint s 
from current data sets JViel et alJ|2006t ISeliak et al.ll2006f) . 

The use of better initial conditions should be a use- 
ful tool in reaching future goals of simulating the nonlinear 
power spectrum to an accuracy of better than 1% percent 
needed f or the next-generation of weak gravitational lensing 
surveys jHuterer fc Takadall2005h . 

The code we used to generate ZA & 2LPT ini- 
tial conditions in this paper is publically available at 
http : // cosmo . nyu . edu/ roman/2LPT 
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APPENDIX A: SECOND-ORDER 
LAGRANGIAN PERTURBATION THEORY 
(2LPT) 

Here we present the basics of 2LPT, we refer the 
reader to the literat ure I^ucherjh99jj [Buch^rt e^aljh 9941 : 
iBouchet et"ailll995l : ICatelarJll995l: IScoccimarrol l2000rj) for 

more information. A detailed step by step implementation 
of 2LPT initial conditions for numerica l simu lations is pre- 
sented in Appendix D2 of IScoccimarrol il998t) . Here we as- 
sume Q m = 1 and S7a = 0. For precision work one must take 
into account the cosmological dependence, including radia- 
tion as well. 

In Lagrangian dynamics, particle positions are de- 
scribed by a displacement field 9 so that, 

x = q+*, (Al) 

where particle trajectories obey the equation of motion, 

g +W (r)£=-V*, (A2) 

where <E> is the gravitational potential. This leads to the 
following equation of motion for the displacement field, 

j2. 



J(q,r)V- 



^ X + H( ) ^ X 
dr 2 dr 



where we have used the Poisson equation and l + <5(x) = 1/J 
where J is the Jacobian of the mapping, 

J (q, t) = Det(5ij (A4) 

and 9 it j = d^i/dqj. 

Solving Eq. I|A30 perturbatively leads to positions and 
velocities, 



A( 2 ) 



-aHV c 



-a HV C 
7 



A(2) 



(A5) 
(A6) 



where the potent ials satisfy the Poisson equations 
jBuchert et alJll994lL 

«5(q) (A7) 

{^(^^(^-[^(q)] 2 }- (as) 



V^ (1) (q) 
V^ (2) (q) 



i>3 



Setting </> (2) = in Equations 1A5IA6I leads to the 
IZel'dovichl lll970D approximation. 

Going beyond 2LPT to third-order (3LPT) becomes 
more costl y due to the need to solve three addi tional Poisson 
equations (iBuchert et alJll994l : rCatelarJll995l) . and the im- 
provement in higher-order statistics is modest (IScoccimarrol 
2000b), leading to an improvement by only a fac tor of two in 
the s cale factor needed to suppress transients (IScoccimarrol 
1998), whereas going from ZA to 2LPT the improvement is 
more than one order of magnitude. However, 3LPT does pro- 
vide a better behavior in underdense regions flBouchet et alJ 
Il995h . and may be worth considering for studies sensitive to 
the statistical properties of voids. 



APPENDIX B: ON GLASS INITIAL 
CONDITIONS 

In this paper, our numerical simulations were started from 
a grid. An alternative to a grid start is to displace parti- 
cles from a configuration known in the literature as a glass, 
which is obtained starting from a Poisson distribution and 
running a n N-body code with the sign of the acceleration 
reversed (ICarlberg fc Couchmanl Il989t iBaueh et alJ Il995l : 



ICouchman et alJll995HWhitdll996USmith et alJl2003 

In this appendix we derive analytically how perturba- 
tions evolve towards a glass, and what steps must be taken 
to avoid exacerbating the problem of transients. Reversing 
accelerations means that the equations of motion are still 
given by Eq. J3J but with 



-1 
3/2 1/2 



(Bl) 



instead of Eq. JIJ. As a result the linear propagator, Eq. JHJ, 
reads 



g a b(v) = e 



/V23 \ 



1 . /V23 \ 

H — ;= sin ri 

v/23 V 4 V 



1 

1 

1 

-6 - 



(B2) 



(A3) 



instead of Eq. © . The "growth" factor D for density pertur- 
bations can be obtained by specifying how initial conditions 
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Figure Bl. Evolution of perturbations towards a glass. The top 
panel shows the growth factor D as a function of scale factor a for 
standard procedure (solid) and extra-damped procedure with u = 
3 (dashed). The bottom panel shows the corresponding behavior 
of the third-moment in terms of its Poisson value. 



are set before evolving towards a glass. Assuming the com- 
mon choice of zero initial velocities, we have D — gu ■ (1, 0) 
and then 



D = 



1 



V23, \ ,1 . /V23, 
— - — In a M = sm In 

4 I \/23 V 4 



(B3) 



where we have switched to the scale factor as a time variable 
(recall rj = In a). Therefore, both eigenmodes contribute an 
overall decaying factor of a _1//4 with oscillations of opposite 
phase. The top panel in Fig. lBll shows Eq. 1B3II as a function 
of a (solid), showing that the decay of perturbations due 
to repulsive gravity is very slow. These damped osci l lation s 
explain the behavior seen in Fig. A2 in lBaueh et alJ lll995l) . 

Using Eq. 11511 with the propagator in Eq. !B2t and 
Gaussian Poisson initial conditions we can derive the be- 
havior of the third moment as evolution proceeds towards a 
glass, 



} glass — 



3D 2 



299^a 
+v / 23(17< 



52 + (23a 



1/4 



75) cos ( 



In a 



,V4 



m ■ / V 23 

9) sm ( - In a 



4, (B4) 



where ap denotes the initial Poisson variance. This result 
is shown as a function of a in the bottom panel in Fig. IB1I 
(solid). Note that (mostly negative) skewness is developed 
as evolution to the glass takes place, therefore although the 
amplitude of perturbations mostly decreases as things evolve 
to a glass, skewness is generated. Unless one evolves by many 
(> 10 2 ) scale factors this can potentially enhance the prob- 
lem of transients. The reason is that if the glass itself has 
negative skewness, this adds to the incorrect skewness that 



initial conditions may have, as we now discuss. This issue 
does not arise for grid initial conditions. 

Particles are displaced from their glass positions q by 
the displacement vector \I> so that their positions are given 
by Eq. 1A1L and the induced density perturbations then 
obey 



(l + 5(q))d 3 g= (l + 5(x))d 3 



(B5) 



where 5(q) are the density perturbations due to the distribu- 
tion of starting positions. These vanish for a grid below the 
Nyquist frequency, but are non-zero for a glass. Then, lin- 
earizing, the total perturbation imposed on a set of particles 
is 

f(x) = (j - l) + 5(q) = 5 ic (x) + <5 g i a ss(x), (B6) 

where the first term is the standard initial condition pertur- 
bation generated by the displacement field, with J given by 
Eq. <A4t , and the second term is the glass perturbation, the 
outcome of the process described by Eqs. 1B2IB41 . Since the 
two perturbations are independent, requiring that the glass 
skewness is negligible (and thus has no impact on transients) 
means to require it to be small compared to the Poisson noise 
of the particles (which typically is smaller than the physical 
perturbations 5i c imposed by the initial condtions). Then 
one must simply require 



(5 )glass <C Up, 



(B7) 



which according to Fig. lBll requires a > 10 2 in the evolution 
towards a glass. 

In the desire for faster convergence towards a glass, the 
Hubble drag term in the equ ations of motion is som etimes 
increased to avoid oscillations llCouchman et alll9 95) . Mak- 
ing such a modification in Eq. 1BH to allow for a generic 
drag of magnitude v gives, 



a, 



-1 
3/2 v 



(B8) 



which leads to propagator eigenmodes time dependence as 
e™ ±?7 with 2 w± — —v ± \/v % — 6, and thus the requirement 
of no oscillations in approaching a glass gives v 2 ^ 6. Tak- 
ing v — 3 as an example, we obtain the growth factor and 
third moment as before, shown in Fig. IB1I as dashed lines. 
We see from this that although the convergence to a glass is 
improved in the growth factor sense, the generation of skew- 
ness is more significant, and therefore one must wait roughly 
the same as with v = 1/2 to have negligible impact. 

Summarizing, generating a glass with an overall expan- 
sion factor by a > 10 2 should be a safe configuration of 
points for starting cosmological simulations. 

This paper has been typeset from a Tj^X/ I^IgX file prepared 
by the author. 
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